function [h_out] = h_func(v,mu,F,h,H,l0,l1,MGF_J)

h0 = v*mu + 0.5*v*h*v'+ l0'*(MGF_J(v)'-1);

h_k = v*F + (MGF_J(v)-1)*l1;

for i = 1:size(H,3)
    h_k(i) = h_k(i) + 0.5*v*H(:,:,i)*v';
end

% h_k(4) = h_k(4) + 0.5*v*H(:,:,4)*v';

h_out = [h0, h_k(2:end)];

end

